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Using Approximating Polynomials in Partial-Global Dynamical 
Simulations 
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Smeared link fermionic actions can be straightforwardly simulated with partial-global updating. The efficiency 
of this simulation is greatly increased if the fermionic matrix is written as a product of several near-identical 
terms. Such a break-up can be achieved using polynomial approximations for the fermionic matrix. In this paper 
we will focus on methods of determining the optimum polynomials. 



1. MOTIVATION 

Dynamical fermion simulations with smeared 
links are used more and more frequently. Smeared 
links improve the chiral and flavor symmetry of 
Wilson and staggered fermionic actions (see ref- 
erences in [Q). The main obstacle in using these 
actions is finding efficient algorithms to simulate 
them. In this paper we present some techniques 
that are used to improve the efficiency of these 
simulations. We will focus here on the HYP ac- 
tion |l|,||,^,D, but we believe that these methods 
are more general and can be used for different 
actions as well. 

The HYP action [| is a modified staggered 
fermion action 



S = S G (U) + S G (V) + S F (V) 



(1) 



where So is the traditional pure gauge action de- 
fined in terms of thin links U, Sq is a gauge action 
defined in terms of smeared links V and 



S F 



-H^lndetO, 



(2) 



where il = M'M restricted to even sites, is the 
staggered fermion action defined in terms of HYP 
links V rather than thin links. The role of So is 
discussed in |J. 

To simulate this action we use a two step ap- 
proach. We first update a part of the thin links, 
U — > U' , using a heath-bath and/or an over- 
relaxation step, then we smear the thin links us- 
ing the hyper-cubic (HYP) blocking || and com- 
pute the fat link action, S G (V) + S F (V). We 



accept this new configuration with probability 
P acc = imn{l, e -^0+^)-M<>+W)} . (3) 

The most difficult step in the algorithm is com- 
puting the determinant ratio 



-S F (v')+S F (v) 



det n (V) 



det ft (V) 
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(4) 



Instead of computing the ratio itself we use a 
stochastic estimator suggested by the last part 
of the relation above 



(5) 
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where A = (O'-^QQ'" 1 / 2 ) This estima- 
tor averages to the determinant ratio since A is 
hermitian and positive definite. Furthermore, the 
variance of this estimator is 



a 2 {E A 
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det (24-1) {det AY 



(6) 



The relation above holds only if the spectrum of A 
is bounded from below by 1/2 otherwise the vari- 
ance diverges. This condition is not necessarily 
fulfilled by our matrix. To address this problem 
we modify our stochastic estimator 



e 'a [f] 



(7) 



This estimator has a finite variance if the low- 
est eigenvalue of A is greater than 1/2™. This 
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condition can be fulfilled if we choose n properly. 

However, this forces us to compute the n th root 

of fermionic matrices. This is why we are forced 

to use the polynomial approximation. 

In order to make the algorithm more efficient 

we employed a reduced matrix [§J^,§|: £l r = 
f 2e -2/(n) j where j 

is a cubic function. The 
fermionic part will be Sp = -f In det O r and the 
reminder is absorbed in Sq- 

Lastly, we mention that it can be proved that 
the algorithm satisfies the detailed balance condi- 
tion for either of the estimators presented above. 

2. POLYNOMIAL APPROXIMATION 

To compute the estimator (Q), where A = 

(tlr 1 ^ 2 Q r Q' r l l 2 ^j , we need to find polyno- 
mial approximations for the matrix functions in- 
volved in the estimator 

P^ip) ~ F(Q)- 1/2n , 

q™(si) ~ F(n) 1/n , (8) 

where F (Q) = (Q exp [-2/ (f2)])™ //4 , n is the 
break-up level and k is the polynomial order. 

For a polynomial to approximate a matrix func- 
tion it needs to approximate the function well 
on the entire spectrum of the matrix. To de- 
termine these polynomials we employed a least 
square method || . The coefficients U of the poly- 
nomial T approximating the function g are deter- 
mined by minimizing the quadratic "distance" 

5 2 (U)= ( 1 dx p(x) (T (x) - g (x) f , (9) 

where Ao, Ai are the spectral bounds of £1 and p 
is a weight function, ideally fi's spectral density. 

This method has the advantage that it turns 
into a linear problem and the solution is 

k 

where k is the polynomial order and 

K iS = f 1 dxp{x)x i+j , 

/>Ai 

9i = / dx p (x) x l g (x) . (11) 

J\ 



For our polynomials we used Ao = 4m 2 , the 
lowest eigenvalue for the staggered fermions ma- 
trix. We set the upper bound to Ai = 16.4 which 
we found to be a good limit for the range of pa- 
rameters we used in our simulations. It should 
be mentioned that the only hard upper bound 
for the spectrum of $1 is 64 + 4m 2 but this limit is 
reduced dynamically. The smearing reduces the 
value of the highest eigenvalue even further. 

The weight function p (x) should approximate 
the spectral density for 17. However, we found 
that for reasonable weight functions the coeffi- 
cients depend very little on the detailed form. 
Thus, we used a weight function of the form 
p (x) = x a which is more convenient for our com- 
putation. The advantage of this form is that we 
can compute the Kij matrix analytically and also 
it generates a recursion relation for the coeffi- 
cients gi. 

To measure the accuracy of the approximation 
we used the distance function (^) where we set 
p(x) = 1. In Fig.|l| we plot the distance between 
[P(")] 2 ™, [Q^] n and F(x) Tl as a function of 
the polynomial order. We notice first that is 
a much better approximation than for n > 1. 
This is due to the fact that F (x) (which is the 
function that approximates) has an almost 
polynomial behavior close to Ao where the func- 
tion is the hardest to approximate. Another inter- 
esting feature that we notice in the graphs is that 
the order of the approximation doesn't improve 
with increasing break-up level. This is puzzling, 
at first, since as n increases the equivalent order 
of the polynomial approximating F (x)^ 1 is n x k. 
However, the number of variables varied to min- 
imize 6 stays the same and it seems that this is 
the determining factor for the order of approxi- 
mation. 

We found that when 5 ~ 10 -8 the numerical 
roundoff error becomes dominant and the higher 
order polynomials do not improve the precision of 
the calculation. Thus, we will be using the poly- 
nomials with S ~ 10 -8 as "exact" ones. We see 
in Fig. |l| that these polynomials have an order 
of k ~ 100. The only exceptions are the polyno- 
mials for n = 1 but they are of no real interest 
since they don't satisfy the minimum eigenvalue 
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Figure 1. The accuracy level for the P and Q 
polynomials: open squares n — 1, filled squares 
n = 4, crosses n — 8 (am — 0.10). 
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Figure 2. Cost vs the polynomial order for m = 
0.04 and n = 4, 8. 



we determine the optimum small ones that min- 
imize the computational cost. The cost is given 
by the number of M'M multiplies and it is deter- 
mined by the order of the polynomials used. We 
found that the order of the optimum small poly- 
nomial is independent of the break-up level (see 
Fig.||) and it increases as we decrease the mass. 
We find that for am = 0.10 we need k — 24, for 
am = 0.06 k = 32 and for am = 0.04 k = 48. 



3. CONCLUSIONS 



requirement. 

To further reduce the cost of computing the 
estimator we implemented a two step approach. 
Since we only use the value of the estimator in the 
accept reject step we only need to know whether 
the value of the estimator is larger or smaller 
that lnr, where r is the random number. We 
will first compute the estimator approximatively, 
using small order polynomials; if the distance be- 
tween the approximate value and lnr is greater 
than e we use this value in the accept reject step. 
If not, we "fall-back", i.e. we compute the "ex- 
act" value of the estimator using the large poly- 
nomials. The value of e is determined in a small 
run as the maximum difference between the ap- 
proximate estimator and its "exact" value. The 
rate of fall-back is usually around 10% and this re- 
duces the computational cost substantially. Hav- 
ing fixed the large polynomials by the S condition 



We have showed how to use the polynomial ap- 
proximation for simulating dynamical fermions. 
First we determine the minimum break-up level n 
using the minimum eigenvalue requirement . Then 
we determine the large polynomials using the 
5 ~ 10 -8 conditions. We then determine the opti- 
mum small polynomials using the minimum cost 
condition. For a more complete discussion see [|J. 
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